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This  paper  is  concerned  with  spectral  estimation  of  a 
finite  number  of  two  dimensional  sinusoids  embedded  in 
white  noise.  Closed  form  expressions  are  derived  for  esti¬ 
mates  using  the  autoregressive  (AR)  prediction  error 
filter  approach,  as  well  as  using  periodogram  with  Bartlett 
window,  and  the  maximum  likelihood  (ML)  method  These 
expressions  are  useful  in  the  study  of  resolving  closely 
spaced  sinusoidal  signals  Over  a  narrow  frequency  band, 
direct  decimation  can  b'5  applied  to  improve  resolution 
and/or  to  reduce  computation.  Simulation  results 
demonstrate  that  decimation  by  ( D\,Dz )  with  a  support 
size  (N\,Nz)  yields  approximately  the  same  resolution  as  a 
support  size  (DXN i,DgN9)  used  with  the  ur.dccimated  sig¬ 
nal  The  use  of  decimation  also  reduces  significantly  com¬ 
putation. 

I.  Introduction 

Computation  rate  and  the  ability  to  resolve  closely 
located  spectral  components  are  of  concern  to  almost  all 
spectral  estimation  methods.  These  problems  have 
received  considerable  attention  in  the  literature  [1-4]  In 
this  paper,  we  extend  some  of  these  results  to  the  two 
dimensional  case.  Specifically,  we  investigate  spectral 
estimates  of  a  finite  number  of  sinusoids  embedded  in 
white  noise  in  two  dimension.  We  shall  concentrate  our 
discussion  on  the  use  of  autoregressive  (AR)  prediction 
error  filter  approach.  Similar  results  can  be  derived  for 
the  maximum  likelihood  (ML)  estimates  and  the  periodo- 
gram  using  Bartlett  window 


II. #  Two-Dimensional  AR  Spectral  Estimation 

The  two-dimensional  process  under  study  is  a  sam¬ 
pled  homogeneous  (stationary)  random  field  |*(ni,nj)J. 
Its  autocorrelation  function  is  defined  as 

r(ni.  na)  =  fijrffc,  +  n,.  fc,  +  ne)x'{kt,  *,)]  (2  1) 

where  E  denotes  expectation  and  •  indicates  complex  con¬ 
jugate  The  power  spectral  density  is 

E  E  r(ni,  n,)  e  *  "**' 

nta’m  n2m  — 

-7TS(,  {STT  (2.2) 

In  practice,  one  observes  far(n i,  ne)l  over  a  finite 
support:  1  s  rij  £  Lt,  1  £  nz  £  /<e  An  estimate  of  the 
autocorrelation  function  can  be  calculated  based  on  the 
observed  data,  and  power  spectral  density  estimate  is 
then  obtained  For  simplicity,  we  shall  use  the  same  nota¬ 
tion  r(nJ(  n*)  and  P((t  ()  to  denote  these  functions  as  well 
as  their  estimates 

Let  (Afj,  N «)  defines  a  rectangular  support  over  which 
the  autocorrelation  function  r(u|.  n8)  is  estimated,  It  is 
convenient  to  define  the  following  notations. 
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C/f  =  [1,  e’\  efZ\  •  .  e  1  j  , 
it  izt  >(A 

C/{  =  [i,  eJ(k  e1*',  ,  e  2  j  , 

V-U<^Vv  (2.3) 

where  ®  denotes  the  direct  product  [5],  and  the  super¬ 
script  T  denotes  transpose.  The  A -A*  co)l mn  vector  L 
can  be  written  as 

V  =  [(/[,  e«Uj,  ,»«[/[.  •  •  .  (2.4) 

That  is,  its  k1*  element,  0  s  k  <  N\Nz,  is  where 

fc~nAg+Z  with  0  5  n  <  N\,  0  i  l  <  Ag.  It  is  con¬ 
venient  to  use  the  two  indices  (n,  l)  rather  than  the  sirglc 

index  k  We  shall  say  that  the  (n,  l)tfi  element  of  V  is 

C/(n  i)  =  0sn<Alt  0s£<A?  (2  5) 

even  though  V  is  a  column  vector. 

Let  X  be  the  column  vector  of  size  A'j.Vg  whose  kiK 
element  is  z(n0+n.  Z0+ £).  where  the  index  pair  (n.  /)  is 
related  to  k  as  b«  We.  and  (n0.  Zo)  are  arbitrary  Wc 
define  the  aulocorrt  atior  matrix  H  as 

R  =  £-[jkA"r]  (2  6) 

Thus  the  element  of  R  ot  the  (n.i)01  row  and  the 
(m ,  k),K  column  is 

=  r(n-m,  l-k). 

0*n ,m <N\.  Q*l,k<Nt  (27) 

Using  the  autoregressive  (AR)  prediction  error  filter 
method  [2],  [6]  the  signal  is  assumed  to  fit  an  AR  model  of 
order  (/V,-l.  Nt~l)  driven  by  a  white  noise  u(n).n*)  It 
can  be  written  as 

JV,-l  Afj-I 

*(n,.nz)  =  -£  E  “u  x(n,-k,  nt-l)  +  u(n,,  ni)2.8) 
t*o  1*0 

where  the  double  summation  does  not  include  the  A:  =1=0 
term  The  coefficients  Oa  are  estimated  by  minimizing 
the  one  step  prediction  error 

A’,-!  A',- 1 

|e(T»!.  n,)|*  =  |*(n,.  n*)  +  E  E  o*t  z(n,-k,  ne-() 2(9) 

Jfc-0  1-0 

This  minimization  leads  to  the  normal  equation 

RA*  =  ep  t  (2  10) 

where  R  is  the  autocorrelation  matrix  given  by  (2.6),  A*  is 
a  AtAg  column  coefficient  vector  whose  fcAg+Z  element  is 
au  with  aoo-l »  €  is  the  A1A2  column  vector 

e  =  [1.  0  0]r  (2  11) 

and  tp  is  the  prediction  error  power,  a  scalar.  The  spec¬ 
trum  is  given  by 

P** «•  O  =  A.-i  V-,  ** -  (2  12) 


E  E 

A- 0  1-0 
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and  can  be  expressed  as 

P*J>«.  0  = 


tTR-'i 

l’rR'c\ 


(2  13) 


where  U  is  given  by  Eq.  (2.3). 

It  is  worth  noting  that  the  2-D  autocorrelation  matrix 
R  defined  above  is  a  symmetric,  positive  definite  snd 
blockToaplxiz ,  but  not  Toeplitz 

Ql.  Sinusoidal  Signals  in  White  Noise 

In  order  to  study  the  resolution  characteristics  of  an 
AR  spectral  estimation,  we  assume  that  true  values  of  the 
autocorrelation  function  ore  known  rather  than  those 
obtained  from  actual  data  The  signal  under  study  is  com¬ 
posed  of  a  finite  number,  K,  sinusoids  and  a  white  noise 
with  unit  power.  The  autocorrelation  function  is 

r(n, .«,)  =  «(»,.  nf)  +  £  W‘)  (3.1) 

*-i 

where  ((*.  (,)  is  the  frequency  of  the  Jfc1*  sinusoid  and  pt 
the  corresponding  power.  The  matrix  R  on  a  support 
(A/ i.  N t)  has  the  form 

/?  =  /  +  £  p,  v;ul  (3.2) 

t-1 

where  /  is  the  (Af ,Ng)  square  identity  matrix  and  U„  is  a 
N,Nt  column  vector  identical  to  U  of  Eq.  (2.3)  but  with 
((».  <*)  in  place  of  ((.  () 

ft  can  be  shown  that  the  AR  spectrum  in  this  case  is 
given  by 

K 

i-En 

P*{<-  ()  =  - - - ^ - (3.3) 

1 1  -  N,N,  £  d,  ?«.*,(<-{«.  f-<i)l* 

where  <U  are  constants  and 
()  = 

with  Bn( A)  given  by 

sin(N\/2) 

N  sin  (A/  2) 

fn  the  case  of  A"=l, 

r  =  /  +  pu;v{. 

and  Eq  (3.3)  reduces  to 

1  -p/  {\*NlNtp) _ 


BhW=  = 

n  n-0 


(3.  ♦) 


(3.5) 


(3.6) 


P«(f.  <)  = 


I  i-fv,/v  <-*,)/  (i+ jv,yv^)  | « ■ 

(3.7) 

It  has  a  peak  at  the  unbiased  location  ({,.  fi),  and  the 
peak  value  is 

PM« i.  «,)  =  (l+AflAT1p)(l+(Af1Af«-l)p)  ±  (Af,Af«p)*(3.8) 

which  is  proportional  to  p*.  So  the  peak  of  AR  spectrum  is 
not  a  power  estimate  but  a  square  power  estimate. 

We  now  determine  the  3  dB  contour  around  the  peak 
in  the  frequency  plane  «,  {)  from  the  equation 

PjU>«,.(,)  =  2  Pm«.0  (3  9) 

By  using  first  2  terms  of  the  Taylor  series  expansion  of  Eq 
(3.3).  it  con  be  shown  that  the  contour  is  approximately 
given  by 

+  {Nt-\){(-(,)\  =  2/  NtNtP  (3  10) 
which  is  a  rhombus  with  the  "major/mir.or  axes”  equal  to 
D(.  ut  =  4/  (N,-\)N,Ngp.  • 

and 

DtAK  =  (3  11) 


This  is  plotted  in  Fig  I.  along  with  some  simulation 
results.  The  data  size  is  L,  =  Lt  =  64  ,  and  the  relevant 
parameters  are:  N,  =  Ng  -  5  and  {i  =  (i  =  C  5n 

For  two  sinusoids  (Af=9'  in  white  noise,  Eq.  (3.3) 
becomes 


PjuiU.t) 


1  —d,—dg 


U-A/.A 'ed.yy ,«„(«, .  IS,)  -  B ,Bgd^n  ^^ag.  (Jt)|* 

(3  12) 

where  a,=(-(,.  and  0i=(-(i.  Pz=(-(g  and  (d,| 

depend  on  the  signal  powers  |pt|  as  well  as  the  frequency 
separations  ({j—  and  (ft— {«).  These  expressions  show 
that  the  PM  are  not  linear  with  respect  to  the  individual 
components,  and  that  there  is  always  interference 
between  them.  The  effect  of  this  interference  on  resolu¬ 
tion  is  not  obvious.  Roughly  speaking,  however,  when  the 
two  frequency  components  are  close  to  each  other  with 
respect  to  the  3  dB  axes,  then  the  two  spectra!  peaks  will 
merge.  Since  the  3  dB  axes  are  decreasing  function  of  the 
signal  power  as  well  as  the  size  of  support,  increasing  sig¬ 
nal  power  and/or  increasing  the  size  of  support  will 
improve  resolution. 

We  note  in  passing  that  if  7Vg=  1  and  f =0,  the  above 
analysis  reduces  to  a  one-dimensional  case  ( 1],  [4]. 


IV.  Other  Spectral  Estimates 

The  above  discussion  can  be  modified  in  a  straight¬ 
forward  way  to  be  applicable  to  the  periodogram  using 
Bartlett  window  PB  [7]  and  the  maximum  likelihood  esti¬ 
mates  Pm  [0].  It  can  be  shown  that  these  estimates  are 
given  respectively  by 


1 


NfN§ 


Qrr* 


Phl«.()  = 


1 


UtR'U’ 


(4.1) 


(42) 


For  K  sinusoids  in  white  noise,  these  expressions 
reduce  to 


and 

Pul«.()  = 


Pa«. «)  =  vV-+  £  (4.3) 

/v  I  /Vg  ft  a  1  *  * 

1  /  NyNt 

K  K 

\-N\Ni^  £  *{am.  0m) 

!•»  mt i  1  e 

(4  4) 

The  coefficients  Q*,  depend  on  the  signal  powers  and  the 
frequency  separations. 


V.  Decimation  to  Improve  Resolution 

It  is  seen  in  the  previous  section  that  the  resolution 
can  be  improved  by  increasing  the  size  of  support.  *n  the 
2-D  case,  however,  the  increasing  size  of  support  will 
greatly  increase  the  computation  since  the  size  of  auto¬ 
correlation  matrix  R  is  {N\Nt)  x  (N\Nz).  We  demonstrate 
in  this  section  that  direct  decimation  of  an  input  data 
sequence  can  improve  the  resolution  without  increasing 
the  size  of  support.  This  technique  has  been  used  in  the 
one-dimensional  case  on  the  periodogram  [3],  the  ML 
(Capon)  method,  and  the  A r  method  [4].  It  is  to  be 
expected  that  the  saving  in  computation  in  the  2-D  case  is 
even  more  significant. 

The  direct  decimation  scheme  is  depicted  in  Fig.  2. 
The  two  reasons  for  the  higher  resolution  are  that  the 
decimation  expands  the  frequency  scales  by  the  factor  D\ 
end  Dt  respectively  and  that  the  bandpass  filter  removes 
interference  from  out-of-band  noise,  thus  raising  the 
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m 
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effective  SNR. 

To  demonstrate  the  improvement  of  resolution  by 
decimation.  Fig  3  shows  the  AR  estimates  on  the  data 
which  consists  of  2  sinusoids  at  frequencies 
(0  1775  r.  C  1775  n)  and  (0  1975  rr.  0  1975  r)  and  noise  7 
dB  below  the  sinusoids  Fig  5(A)  shows  the  P ax  with 
A^i  =  =  5  without  decimation,  (B)  the  Par.d  with 
Nx  =  Ni  =  5  and  decimation  D\  =  Dt  =  4,  both  plotted  on 
0  126  r  <  s  025rr.  These  simulation  results  show 
that  direct  decimation  improves  the  resolution  without 
increasing  the  si2e  of  support 

Ir.  order  to  investigate  the  improvement  of  resolu¬ 
tion  by  decimation,  we  may  introduce  the  notation  of 
resolution  boundary"  which,  as  in  the  1-D  case  [9],  for 
twn  •  .usoids  with  equal  power  is  defined  as  the  frequency 
separation  (A<.  A<)  =  Ui-$»l)  at  which  the  spec¬ 

trum  at  the  center  frequency  is  equal  to  the  average  of 
the  spectra  at  the  two  sinusoid  frequencies,  i.e., 

=  | (i)  +  ?*)]  (5.1) 

The  resolution  boundary  is  the  minimum  resolvable  fre¬ 
quency  separation  for  a  given  SNR.  In  Fig  4  we  show  a 
special  symmetric  case  with  The  solid  line  exhi¬ 

bits  the  theoretical  curve  of  SNR  vs.  computed  from 
(5  1)  with  N\  =  =  io  The  dashed  line  indicates  the 

same  theoretical  curve  for  corresponding  decimated 
spectra  with  Nx  -  A’z  =  5,  and  Dx  =  Dz  =  2.  The  two  curves 
are  seen  to  be  close  to  each  other  The  triangles  and  rec¬ 
tangles  are  corresponding  simulation  results  on  a 
(32  x  32)  data  with  center  frequency  =  0.25  rr. 

Since  these  results  are  obtained  by  averaging  only  two 
independent  trials,  a  considerable  variation  is  present. 


VL  Decimation  to  Reduce  Computation 

In  the  last  section,  we  have  seen  that  using  decima¬ 
tion  by  factor  (Dlt  Dz)  can  reduce  the  support  size 

(N\,  Nz)  to  (-7^-.  ~-)  while  maintaining  the  same  resolu- 
U\  U  a 

tier..  A  reduction  in  the  size  of  support  is  accompanied  by 
a  saving  of  computation. 

For  a  support  (A,,  Afg)  and  a  data  size  {Lx  x  Lz),  the 
number  of  multiplication  in  the  computation  of  autocorre¬ 
lation  matrix  is 

*  N,N,L,L,  (6.1) 

If  the  Justic  algorithm  [10]  is  used  to  solve  the  block 
Toeplitz  matrix  normal  equation,  Eq.  (2.10),  the  computa¬ 
tion  may  be  0 (A*A|).  If  the  Gaussian-Seidel  iteration  is 
used,  the  number  of  multiplication  is 

Ms  *  ?(A \NZ)Z  (6.2) 

where  7  is  a  constant  depended  on  the  number  of  itera¬ 
tion  Thus  the  total  number  of  multiplication  for  a  Pas 
estimate  is 

««  =  Ur  +  Un  «*  N,N,L,L,  +  jNfNt  (6.3) 

Suppose  the  band-pass  filter  preceding  the  down 
sampling  in  Fig  2  is  a  first  quadrant  FIR  filter  with  length 
of  impulse  response  (Ay,,  Ay,)  The  number  of  multipli¬ 
cation  for  the  filtering  is 

..  _  (ti+/'V,)AVi  (Lj+AVfJAy, 

Urg - ZD, - ZD, -  <64) 

The  number  of  multiplication  for  computing  the 
autocorrelation  matrix  and  for  solving  normal  equation 
for  a  decimated  AR  spectrum  are  respectively 

u  —  Lj+Nn  N i+Di  (  Lg+Nrt  A/g+Z?e  ] 

"K  0~  D,D,  D,  ~  ZD,  Dt  ZD, 


(Af.Af,)*  U„ 

ard  Mk  d~7  (W^  {66) 

Thus,  the  total  number  of  multiplication  is  approximately 
Max  d  =  Mr.  d  +  Ms ,d  +  Mf.D  =  \z~+  ^f.d  (6.7) 

\U\UZ, 

If  D\  -  Dt  -  D,  the  saving  of  computation  is  by  a  factor 
D*.  excluding  the  filtering  Mf.d 


Discussion 

The  resolution  of  sinusoids  in  white  noise  using  2D  AR 
spectrum  is  investigated  in  this  paper.  Closed  form 
expressions  of  spectral  estimates  are  given.  The  peaks  of 
the  AR  spectra  indicate  the  square  of  the  power  of  each 
component.  A  3  dB  contour  in  the  frequency  plane  is 
introduced  to  facilitate  the  study  of  resolution  charac¬ 
teristics. 

Decimation  is  then  applied  to  narrow-band  2-D  spec¬ 
tral  estimation.  Simulation  results  indicate  that  a  spec¬ 
tral  estimate  on  a  support  (NXNZ)  with  a  decimation  fac¬ 
tor  of  (Z?|,  Dt)  has  a  resolution  approximately  equivalent 
to  that  of  a  spectral  estimate  on  a  support  (DXNX,  DzNz) 
without  decimation.  It  is  also  shown  that  decimation 
reduces  computation  by  a  factor  (D\DZ)Z  without 
sacrificing  resolution 

This  analysis  can  be  extended  to  most  other  2-D 
spectral  estimation  methods,  such  as  periodogram  and 
maximum  likelihood  method,  and  similar  results  car  be 
obtained 
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